/**
 * @file   radixsort.h
 * @author karpar Hu <huzhili@gmail.com>
 * @date   Fri May  1 21:38:36 2009
 * 
 * @brief  radixsort routines,for description of this see "Enginerring Radix Sort"
 *         by D.McIlroy and wikipedia link:http://en.wikipedia.org/wiki/Radix_sort
 *         which is more in detail.It faster than quicksort by 3 times, when soring 
 *         string.
 */

/**
 * An example
 * Original, unsorted list:
 *   170, 45, 75, 90, 2, 24, 802, 66
 * Sorting by least significant digit (1s place) gives:
 *   170, 90, 2, 802, 24, 45, 75, 66
 * Notice that we keep 2 before 802, because 2 occurred before 802 in the original list,
 * and similarly for pairs 170 & 90 and 45 & 75.
 * Sorting by next digit (10s place) gives:
 *   2, 802, 24, 45, 66, 170, 75, 90
 * Sorting by most significant digit (100s place) gives:
 *   2, 24, 45, 66, 75, 90, 170, 802
 * Most significant digit radix sorts
 * Recursion
 * A recursively subdividing MSD radix sort algorithm works as follows:
 * 1.Take the most significant digit of each key.
 * 2.Sort the list of elements based on that digit, grouping elements with
 * the same digit into one bucket.
 * 3.Recursively sort each bucket, starting with the next digit to the right.
 * 4.Concatenate the buckets together in order.
 * 
 */

#ifndef _RADIXSORT_H_
#define _RADIXSORT_H_

//#include "zdefs.h"
#include <errno.h>

typedef unsigned int u_int;
typedef unsigned char u_char;

typedef struct stack {
  const u_char **sa;            /**< the string array pointer */
  int sn, si;                   /**< sn:the length of each subarray, si:string index */
} stack;

#define empty(s) (s >= sp)
#define pop(a, n, i)  a = (--sp)->sa, n = sp->sn, i = sp->si
#define push(a, n, i) sp->sa = a, sp->sn = n, (sp++)->si = i
#define swap(a, b, t) t = a, a = b, b = t

//static inline void
//simplesort(const u_char **a, int n, int i,  const u_char *tr, u_int endch) __attribute__((always_inline));

static inline void
simplesort(const u_char **a, int n, int i,  const u_char *tr, u_int endch)
{
  u_char ch;
  const u_char **ak, **ai, *s, *t;
  
  for (ak = a + 1; --n >= 1; ak++) 
    for (ai = ak; ai > a; ai--) {
 
     for (s = ai[0] + i, t = ai[-1] + i; (ch = tr[*s]) != endch; s++, t++)
        if (ch != tr[*t])
          break;
      if (ch >=  tr[*t])
        break;
      swap(ai[0], ai[-1], s);
    }
}

static void r_sort_a(const u_char **, int, int, const u_char *, u_int);
static void r_sort_b(const u_char **, const u_char **, int , int, const u_char *, u_int);

#define THRESHOLD 20            /**< divert to simplesort() */
#define SIZE 512                /**< default stack size */

  
#define SETUP {                                 \
    if (tab == NULL) {                          \
      tr = tr0;                                 \
      for (c = 0; c < endch; c++)               \
        tr0[c] = c + 1;                         \
      tr0[c] = 0;                               \
      for (c++; c < 256; c++)                   \
        tr0[c] = c;                             \
      endch = 0;                                \
    } else {                                    \
      endch = tab[endch];                       \
      tr = tab;                                 \
      if (endch != 0 && endch != 255) {         \
        errno = EINVAL;                         \
      }                                         \
    }                                           \
  }                                             \

/** 
 * unstable sort,but in-place.
 * 
 * @param a the string array
 * @param tab charater table
 * @param n the number of elements
 * @param endch the terminate character.
 * 
 * @return 0:all right, 1: error happened
 */
int radixsort(const u_char **a, const u_char *tab, int n, u_int endch)
{
  const u_char *tr;
  int c;
  u_char tr0[256];
  
  if (tab == NULL) {                         
    tr = tr0;                               
    for (c = 0; c < endch; c++)             
      tr0[c] = c + 1;                       
    tr0[c] = 0;                             
    for (c++; c < 256; c++)                 
      tr0[c] = c;                           
    endch = 0;                              
  } else {                                  
    endch = tab[endch];                     
    tr = tab;                               
    if (endch != 0 && endch != 255) {       
      errno = EINVAL;                       
    }                                       
  }                                         

  r_sort_a(a, n, 0, tr, endch);
  return 0;
}


static void
r_sort_a(const u_char **a, int n, int i, const u_char *tr, u_int endch)
{
  static int count[256], nc, bmin;
  int c;
  const u_char **ak, *r;
  stack s[SIZE], *sp, *sp0, *sp1, tmp;
  int *cp, bigc;
  const u_char **an, *t, **aj, **top[256];
  
  /* Set up stack. */
  sp = s;
  push(a, n, i);
  while (!empty(s)) {
    pop(a, n, i);
    if (n < THRESHOLD) {
      simplesort(a, n, i, tr, endch);
      continue;
    }
    an = a + n;
    
    /* Make character histogram. */
    if (nc == 0) {
      bmin = 255;
      for (ak = a; ak < an;) {
        c = tr[(*ak++)[i]];
        if (++count[c] == 1 && c != endch) {
          if (c < bmin) 
            bmin = c;
          nc++;
        }
      }
      if (sp + nc < s + SIZE) {  /* get more stack */
        r_sort_a(a, n, i, tr, endch);
        continue;
      }
    }
    /* Special case: if all strings have the same
     * character at position i, move on to the next character.
     */
    if (nc == 1 && count[bmin] == n) {
      push(a, n , i + 1);
      nc = count[bmin] = 0;
      continue;
    }
    /* Set top[]; push incompletely sorted bins onto stack. top[] = pointers to last 
     * out-of-place element in bins. count[] = counts of elements in bins. 
     * Before permuting: top[c-1] + count[c] = top[c]; 
     * during deal: top[c] counts down to top[c-1]. 
     */
    sp0 = sp1 = sp;             /*  sp1:Stack position of biggest bin. */
    bigc = 2;                   /* Size of biggest bin. */
    if (endch == 0)             /* Special case: set top[eos]. */
      top[0] = ak = a + count[0];
    else {
      ak = a;
      top[255] = an;
    }
    for (cp = count + bmin; nc > 0; cp++) {
      while (*cp == 0)          /* Find next non-empty pile. */
        cp++;
      if (*cp > 1) {
        if (*cp > bigc) {
          bigc = *cp;
          sp1 = sp;
        }
        push(ak, *cp, i+1);
      }
      top[cp - count] = ak += *cp;
      nc--;
    }
    swap(*sp0, *sp1, tmp);      /* put the biggest to last. */
    /* permute in-place, go for papers above to get more detail */
    for (aj = a; aj < an; *aj = r, aj += count[c], count[c] = 0) 
      for (r = *aj; aj < (ak = --top[c = tr[r[i]]]);)
        swap(*ak, r, t);
  }
}



#endif /* _RADIXSORT_H_ */
